library(tidyverse)
library(bayesrules)
library(bayesplot)
library(rstan)
library(rstanarm)
library(grid)
library(gridExtra)
library(janitor)
library(broom.mixed)
library(tidybayes)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ── ✔ dplyr 1.1.2 ✔ readr 2.1.4 ✔ forcats 1.0.0 ✔ stringr 1.5.0 ✔ ggplot2 3.4.2 ✔ tibble 3.2.1 ✔ lubridate 1.9.2 ✔ tidyr 1.3.0 ✔ purrr 1.0.1 ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ── ✖ dplyr::filter() masks stats::filter() ✖ dplyr::lag() masks stats::lag() ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors This is bayesplot version 1.10.0 - Online documentation and vignettes at mc-stan.org/bayesplot - bayesplot theme set to bayesplot::theme_default() * Does _not_ affect other ggplot2 plots * See ?bayesplot_theme_set for details on theme setting Loading required package: StanHeaders rstan (Version 2.21.8, GitRev: 2e1f913d3ca3) For execution on a local, multicore CPU with excess RAM we recommend calling options(mc.cores = parallel::detectCores()). To avoid recompilation of unchanged Stan programs, we recommend calling rstan_options(auto_write = TRUE) Attaching package: ‘rstan’ The following object is masked from ‘package:tidyr’: extract Loading required package: Rcpp This is rstanarm version 2.21.4 - See https://mc-stan.org/rstanarm/articles/priors for changes to default priors! - Default priors may change, so it's safest to specify priors, even if equivalent to the defaults. - For execution on a local, multicore CPU with excess RAM we recommend calling options(mc.cores = parallel::detectCores()) Attaching package: ‘rstanarm’ The following object is masked from ‘package:rstan’: loo Attaching package: ‘gridExtra’ The following object is masked from ‘package:dplyr’: combine Attaching package: ‘janitor’ The following objects are masked from ‘package:stats’: chisq.test, fisher.test
options(repr.plot.width=15, repr.plot.height=4)
plot_normal(mean = 0.08, sd=0.1)
plot_normal(mean = 0, sd=0.02)
x <- seq(0,0.3,0.001)
ggplot( data.frame(x=x, p=dexp( x, rate=6.25 ) ) ) + geom_line( aes(x=x, y=p))
options(repr.plot.width=15, repr.plot.height=4)
plot_normal(mean = 0, sd=5)
plot_normal(mean = 1, sd=0.1)
plot_normal(mean = 5000, sd = 2000)
plot_normal(mean = 10, sd = 5)
bike_humidity_model_prior <- stan_glm(
rides ~ humidity, data = bikes,
family = gaussian,
prior_intercept = normal(5000, 2000),
prior = normal(10, 5),
prior_aux = exponential(1/2000),
chains = 5, iter = 4000*2, seed = 84735, prior_PD = TRUE)
SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1). Chain 1: Chain 1: Gradient evaluation took 1.4e-05 seconds Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds. Chain 1: Adjust your expectations accordingly! Chain 1: Chain 1: Chain 1: Iteration: 1 / 8000 [ 0%] (Warmup) Chain 1: Iteration: 800 / 8000 [ 10%] (Warmup) Chain 1: Iteration: 1600 / 8000 [ 20%] (Warmup) Chain 1: Iteration: 2400 / 8000 [ 30%] (Warmup) Chain 1: Iteration: 3200 / 8000 [ 40%] (Warmup) Chain 1: Iteration: 4000 / 8000 [ 50%] (Warmup) Chain 1: Iteration: 4001 / 8000 [ 50%] (Sampling) Chain 1: Iteration: 4800 / 8000 [ 60%] (Sampling) Chain 1: Iteration: 5600 / 8000 [ 70%] (Sampling) Chain 1: Iteration: 6400 / 8000 [ 80%] (Sampling) Chain 1: Iteration: 7200 / 8000 [ 90%] (Sampling) Chain 1: Iteration: 8000 / 8000 [100%] (Sampling) Chain 1: Chain 1: Elapsed Time: 0.24437 seconds (Warm-up) Chain 1: 0.040126 seconds (Sampling) Chain 1: 0.284496 seconds (Total) Chain 1: SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2). Chain 2: Chain 2: Gradient evaluation took 7e-06 seconds Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds. Chain 2: Adjust your expectations accordingly! Chain 2: Chain 2: Chain 2: Iteration: 1 / 8000 [ 0%] (Warmup) Chain 2: Iteration: 800 / 8000 [ 10%] (Warmup) Chain 2: Iteration: 1600 / 8000 [ 20%] (Warmup) Chain 2: Iteration: 2400 / 8000 [ 30%] (Warmup) Chain 2: Iteration: 3200 / 8000 [ 40%] (Warmup) Chain 2: Iteration: 4000 / 8000 [ 50%] (Warmup) Chain 2: Iteration: 4001 / 8000 [ 50%] (Sampling) Chain 2: Iteration: 4800 / 8000 [ 60%] (Sampling) Chain 2: Iteration: 5600 / 8000 [ 70%] (Sampling) Chain 2: Iteration: 6400 / 8000 [ 80%] (Sampling) Chain 2: Iteration: 7200 / 8000 [ 90%] (Sampling) Chain 2: Iteration: 8000 / 8000 [100%] (Sampling) Chain 2: Chain 2: Elapsed Time: 0.166822 seconds (Warm-up) Chain 2: 0.032306 seconds (Sampling) Chain 2: 0.199128 seconds (Total) Chain 2: SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3). Chain 3: Chain 3: Gradient evaluation took 4e-06 seconds Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds. Chain 3: Adjust your expectations accordingly! Chain 3: Chain 3: Chain 3: Iteration: 1 / 8000 [ 0%] (Warmup) Chain 3: Iteration: 800 / 8000 [ 10%] (Warmup) Chain 3: Iteration: 1600 / 8000 [ 20%] (Warmup) Chain 3: Iteration: 2400 / 8000 [ 30%] (Warmup) Chain 3: Iteration: 3200 / 8000 [ 40%] (Warmup) Chain 3: Iteration: 4000 / 8000 [ 50%] (Warmup) Chain 3: Iteration: 4001 / 8000 [ 50%] (Sampling) Chain 3: Iteration: 4800 / 8000 [ 60%] (Sampling) Chain 3: Iteration: 5600 / 8000 [ 70%] (Sampling) Chain 3: Iteration: 6400 / 8000 [ 80%] (Sampling) Chain 3: Iteration: 7200 / 8000 [ 90%] (Sampling) Chain 3: Iteration: 8000 / 8000 [100%] (Sampling) Chain 3: Chain 3: Elapsed Time: 0.187401 seconds (Warm-up) Chain 3: 0.03235 seconds (Sampling) Chain 3: 0.219751 seconds (Total) Chain 3: SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4). Chain 4: Chain 4: Gradient evaluation took 5e-06 seconds Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds. Chain 4: Adjust your expectations accordingly! Chain 4: Chain 4: Chain 4: Iteration: 1 / 8000 [ 0%] (Warmup) Chain 4: Iteration: 800 / 8000 [ 10%] (Warmup) Chain 4: Iteration: 1600 / 8000 [ 20%] (Warmup) Chain 4: Iteration: 2400 / 8000 [ 30%] (Warmup) Chain 4: Iteration: 3200 / 8000 [ 40%] (Warmup) Chain 4: Iteration: 4000 / 8000 [ 50%] (Warmup) Chain 4: Iteration: 4001 / 8000 [ 50%] (Sampling) Chain 4: Iteration: 4800 / 8000 [ 60%] (Sampling) Chain 4: Iteration: 5600 / 8000 [ 70%] (Sampling) Chain 4: Iteration: 6400 / 8000 [ 80%] (Sampling) Chain 4: Iteration: 7200 / 8000 [ 90%] (Sampling) Chain 4: Iteration: 8000 / 8000 [100%] (Sampling) Chain 4: Chain 4: Elapsed Time: 0.19712 seconds (Warm-up) Chain 4: 0.032916 seconds (Sampling) Chain 4: 0.230036 seconds (Total) Chain 4: SAMPLING FOR MODEL 'continuous' NOW (CHAIN 5). Chain 5: Chain 5: Gradient evaluation took 6e-06 seconds Chain 5: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds. Chain 5: Adjust your expectations accordingly! Chain 5: Chain 5: Chain 5: Iteration: 1 / 8000 [ 0%] (Warmup) Chain 5: Iteration: 800 / 8000 [ 10%] (Warmup) Chain 5: Iteration: 1600 / 8000 [ 20%] (Warmup) Chain 5: Iteration: 2400 / 8000 [ 30%] (Warmup) Chain 5: Iteration: 3200 / 8000 [ 40%] (Warmup) Chain 5: Iteration: 4000 / 8000 [ 50%] (Warmup) Chain 5: Iteration: 4001 / 8000 [ 50%] (Sampling) Chain 5: Iteration: 4800 / 8000 [ 60%] (Sampling) Chain 5: Iteration: 5600 / 8000 [ 70%] (Sampling) Chain 5: Iteration: 6400 / 8000 [ 80%] (Sampling) Chain 5: Iteration: 7200 / 8000 [ 90%] (Sampling) Chain 5: Iteration: 8000 / 8000 [100%] (Sampling) Chain 5: Chain 5: Elapsed Time: 0.312349 seconds (Warm-up) Chain 5: 0.033881 seconds (Sampling) Chain 5: 0.34623 seconds (Total) Chain 5:
library(tidybayes)
options(repr.plot.width=15, repr.plot.height=10)
bikes %>%
add_fitted_draws(bike_humidity_model_prior, n = 100) %>%
ggplot(aes(x = humidity, y = rides)) +
geom_line(aes(y = .value, group = .draw), alpha = 0.1, col = "blue") +
geom_point(data = bikes, size = 0.05)
Warning message: “`fitted_draws` and `add_fitted_draws` are deprecated as their names were confusing. - Use [add_]epred_draws() to get the expectation of the posterior predictive. - Use [add_]linpred_draws() to get the distribution of the linear predictor. - For example, you used [add_]fitted_draws(..., scale = "response"), which means you most likely want [add_]epred_draws(...). NOTE: When updating to the new functions, note that the `model` parameter is now named `object` and the `n` parameter is now named `ndraws`.”
bikes %>%
add_predicted_draws(bike_humidity_model_prior, n = 4) %>%
ggplot(aes(x = humidity, y = rides)) +
geom_point(aes(y = .prediction, group = .draw), size = 0.2) +
facet_wrap(~ .draw)
Warning message:
“
In add_predicted_draws(): The `n` argument is a deprecated alias for `ndraws`.
Use the `ndraws` argument instead.
See help("tidybayes-deprecated").
”
Four samples from exponential distribution for prior on sigma:
rexp(4, rate=1/2000)
ggplot( bikes, aes(x=humidity, y=rides) ) + geom_point() + geom_smooth( method="lm" )
`geom_smooth()` using formula = 'y ~ x'
options(repr.plot.width=15, repr.plot.height=8)
ggplot( bikes %>% filter( humidity > 70, humidity < 80), aes(x=rides) ) + geom_histogram()
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
options(repr.plot.width=15, repr.plot.height=8)
ggplot( bikes %>% filter( humidity > 50, humidity < 60), aes(x=rides) ) + geom_histogram()
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
The variance due to noise is extremely large, so the choice of a low-variance linear model is the best. The distributions of Y along different regions of X look to some approximation normal, so one can try with a normal model, even though this assumption is not perfect. The variance appears to be more or less constant along the x-axis.
bike_humidity_model <- stan_glm(
rides ~ humidity, data = bikes,
family = gaussian,
prior_intercept = normal(5000, 2000),
prior = normal(10, 5),
prior_aux = exponential(1/2000),
chains = 5, iter = 4000*2, seed = 84735, prior_PD = FALSE)
SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1). Chain 1: Chain 1: Gradient evaluation took 1.3e-05 seconds Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.13 seconds. Chain 1: Adjust your expectations accordingly! Chain 1: Chain 1: Chain 1: Iteration: 1 / 8000 [ 0%] (Warmup) Chain 1: Iteration: 800 / 8000 [ 10%] (Warmup) Chain 1: Iteration: 1600 / 8000 [ 20%] (Warmup) Chain 1: Iteration: 2400 / 8000 [ 30%] (Warmup) Chain 1: Iteration: 3200 / 8000 [ 40%] (Warmup) Chain 1: Iteration: 4000 / 8000 [ 50%] (Warmup) Chain 1: Iteration: 4001 / 8000 [ 50%] (Sampling) Chain 1: Iteration: 4800 / 8000 [ 60%] (Sampling) Chain 1: Iteration: 5600 / 8000 [ 70%] (Sampling) Chain 1: Iteration: 6400 / 8000 [ 80%] (Sampling) Chain 1: Iteration: 7200 / 8000 [ 90%] (Sampling) Chain 1: Iteration: 8000 / 8000 [100%] (Sampling) Chain 1: Chain 1: Elapsed Time: 0.136099 seconds (Warm-up) Chain 1: 0.139633 seconds (Sampling) Chain 1: 0.275732 seconds (Total) Chain 1: SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2). Chain 2: Chain 2: Gradient evaluation took 6e-06 seconds Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds. Chain 2: Adjust your expectations accordingly! Chain 2: Chain 2: Chain 2: Iteration: 1 / 8000 [ 0%] (Warmup) Chain 2: Iteration: 800 / 8000 [ 10%] (Warmup) Chain 2: Iteration: 1600 / 8000 [ 20%] (Warmup) Chain 2: Iteration: 2400 / 8000 [ 30%] (Warmup) Chain 2: Iteration: 3200 / 8000 [ 40%] (Warmup) Chain 2: Iteration: 4000 / 8000 [ 50%] (Warmup) Chain 2: Iteration: 4001 / 8000 [ 50%] (Sampling) Chain 2: Iteration: 4800 / 8000 [ 60%] (Sampling) Chain 2: Iteration: 5600 / 8000 [ 70%] (Sampling) Chain 2: Iteration: 6400 / 8000 [ 80%] (Sampling) Chain 2: Iteration: 7200 / 8000 [ 90%] (Sampling) Chain 2: Iteration: 8000 / 8000 [100%] (Sampling) Chain 2: Chain 2: Elapsed Time: 0.094059 seconds (Warm-up) Chain 2: 0.136463 seconds (Sampling) Chain 2: 0.230522 seconds (Total) Chain 2: SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3). Chain 3: Chain 3: Gradient evaluation took 7e-06 seconds Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds. Chain 3: Adjust your expectations accordingly! Chain 3: Chain 3: Chain 3: Iteration: 1 / 8000 [ 0%] (Warmup) Chain 3: Iteration: 800 / 8000 [ 10%] (Warmup) Chain 3: Iteration: 1600 / 8000 [ 20%] (Warmup) Chain 3: Iteration: 2400 / 8000 [ 30%] (Warmup) Chain 3: Iteration: 3200 / 8000 [ 40%] (Warmup) Chain 3: Iteration: 4000 / 8000 [ 50%] (Warmup) Chain 3: Iteration: 4001 / 8000 [ 50%] (Sampling) Chain 3: Iteration: 4800 / 8000 [ 60%] (Sampling) Chain 3: Iteration: 5600 / 8000 [ 70%] (Sampling) Chain 3: Iteration: 6400 / 8000 [ 80%] (Sampling) Chain 3: Iteration: 7200 / 8000 [ 90%] (Sampling) Chain 3: Iteration: 8000 / 8000 [100%] (Sampling) Chain 3: Chain 3: Elapsed Time: 0.094782 seconds (Warm-up) Chain 3: 0.140161 seconds (Sampling) Chain 3: 0.234943 seconds (Total) Chain 3: SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4). Chain 4: Chain 4: Gradient evaluation took 5e-06 seconds Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds. Chain 4: Adjust your expectations accordingly! Chain 4: Chain 4: Chain 4: Iteration: 1 / 8000 [ 0%] (Warmup) Chain 4: Iteration: 800 / 8000 [ 10%] (Warmup) Chain 4: Iteration: 1600 / 8000 [ 20%] (Warmup) Chain 4: Iteration: 2400 / 8000 [ 30%] (Warmup) Chain 4: Iteration: 3200 / 8000 [ 40%] (Warmup) Chain 4: Iteration: 4000 / 8000 [ 50%] (Warmup) Chain 4: Iteration: 4001 / 8000 [ 50%] (Sampling) Chain 4: Iteration: 4800 / 8000 [ 60%] (Sampling) Chain 4: Iteration: 5600 / 8000 [ 70%] (Sampling) Chain 4: Iteration: 6400 / 8000 [ 80%] (Sampling) Chain 4: Iteration: 7200 / 8000 [ 90%] (Sampling) Chain 4: Iteration: 8000 / 8000 [100%] (Sampling) Chain 4: Chain 4: Elapsed Time: 0.152082 seconds (Warm-up) Chain 4: 0.140769 seconds (Sampling) Chain 4: 0.292851 seconds (Total) Chain 4: SAMPLING FOR MODEL 'continuous' NOW (CHAIN 5). Chain 5: Chain 5: Gradient evaluation took 6e-06 seconds Chain 5: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds. Chain 5: Adjust your expectations accordingly! Chain 5: Chain 5: Chain 5: Iteration: 1 / 8000 [ 0%] (Warmup) Chain 5: Iteration: 800 / 8000 [ 10%] (Warmup) Chain 5: Iteration: 1600 / 8000 [ 20%] (Warmup) Chain 5: Iteration: 2400 / 8000 [ 30%] (Warmup) Chain 5: Iteration: 3200 / 8000 [ 40%] (Warmup) Chain 5: Iteration: 4000 / 8000 [ 50%] (Warmup) Chain 5: Iteration: 4001 / 8000 [ 50%] (Sampling) Chain 5: Iteration: 4800 / 8000 [ 60%] (Sampling) Chain 5: Iteration: 5600 / 8000 [ 70%] (Sampling) Chain 5: Iteration: 6400 / 8000 [ 80%] (Sampling) Chain 5: Iteration: 7200 / 8000 [ 90%] (Sampling) Chain 5: Iteration: 8000 / 8000 [100%] (Sampling) Chain 5: Chain 5: Elapsed Time: 0.173074 seconds (Warm-up) Chain 5: 0.142273 seconds (Sampling) Chain 5: 0.315347 seconds (Total) Chain 5:
Trace plot:
options(repr.plot.width=15, repr.plot.height=4)
mcmc_trace( bike_humidity_model )
It looks like the chains have mixed.
Densities:
mcmc_dens_overlay( bike_humidity_model )
Same here. Then probably Rhat is close to 1:
rhat( bike_humidity_model )
Autocorrelation:
mcmc_acf( bike_humidity_model )
Warning message: “The `facets` argument of `facet_grid()` is deprecated as of ggplot2 2.2.0. ℹ Please use the `rows` argument instead. ℹ The deprecated feature was likely used in the bayesplot package. Please report the issue at <https://github.com/stan-dev/bayesplot/issues/>.”
Looks good! Then probably also effective sample size ratio is large.
neff_ratio( bike_humidity_model )
Posterior lines:
options(repr.plot.width=15, repr.plot.height=10)
bikes %>%
add_fitted_draws(bike_humidity_model, n = 100) %>%
ggplot(aes(x = humidity, y = rides)) +
geom_line(aes(y = .value, group = .draw), alpha = 0.1, col = "blue") +
geom_point(data = bikes, size = 0.05)
Warning message: “`fitted_draws` and `add_fitted_draws` are deprecated as their names were confusing. - Use [add_]epred_draws() to get the expectation of the posterior predictive. - Use [add_]linpred_draws() to get the distribution of the linear predictor. - For example, you used [add_]fitted_draws(..., scale = "response"), which means you most likely want [add_]epred_draws(...). NOTE: When updating to the new functions, note that the `model` parameter is now named `object` and the `n` parameter is now named `ndraws`.”
The intercept was narrowed down quite a bit!
Posterior distributions:
data.frame( bike_humidity_model ) %>%
summarize_all( list(mean=mean, sd=sd) )
| X.Intercept._mean | humidity_mean | sigma_mean | X.Intercept._sd | humidity_sd | sigma_sd |
|---|---|---|---|---|---|
| <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> |
| 3433.796 | 0.8106482 | 1579.088 | 226.2688 | 3.403878 | 50.14182 |
The variances have decreased by a large factor! (2000 $\to$ 226, 10 $\to$ 3.4, 2000 $\to$ 50)
tidy(bike_humidity_model, effects = c("fixed", "aux"),
conf.int = TRUE, conf.level = 0.95)
| term | estimate | std.error | conf.low | conf.high |
|---|---|---|---|---|
| <chr> | <dbl> | <dbl> | <dbl> | <dbl> |
| (Intercept) | 3434.7221760 | 225.250742 | 2990.848174 | 3876.236774 |
| humidity | 0.8099946 | 3.422603 | -5.816559 | 7.489818 |
| sigma | 1578.2656423 | 50.141823 | 1483.217699 | 1680.801519 |
| mean_PPD | 3484.2562015 | 99.792313 | 3289.073650 | 3679.794383 |
The noise is a bit lower than expected, however still extremely large.
Plot of prior (red) and posterior distribution (blue) of slope:
options(repr.plot.width=15, repr.plot.height=4)
ggplot() +
geom_histogram( data=data.frame(bike_humidity_model), aes(x=humidity), alpha=0.5, fill="blue", ) +
geom_histogram( data=data.frame(bike_humidity_model_prior), aes(x=humidity), alpha=0.5, fill="red")
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`. `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
We strongly assumed a positive slope in the prior, however the data suggest that the slope could with equal probability be negative and is (given data + prior) most likely near zero.
No, we do not have ample evidence. The probability that the slope is positive is still larger than 50%.
bike_humidity_model_df = data.frame( bike_humidity_model )
head( bike_humidity_model_df )
| X.Intercept. | humidity | sigma | |
|---|---|---|---|
| <dbl> | <dbl> | <dbl> | |
| 1 | 3113.506 | 4.8994985 | 1634.080 |
| 2 | 3260.669 | 2.1678773 | 1595.147 |
| 3 | 3416.404 | 2.6359108 | 1559.959 |
| 4 | 3407.504 | 2.2857136 | 1570.255 |
| 5 | 3489.423 | -1.1614780 | 1582.603 |
| 6 | 3303.120 | 0.8239054 | 1648.416 |
pred90 <- bike_humidity_model_df %>%
mutate( mu = X.Intercept. + humidity * 90 ) %>%
mutate( Ynew = rnorm( nrow(bike_humidity_model_df), mean=mu, sd=sigma ) )
head( pred90 )
| X.Intercept. | humidity | sigma | mu | Ynew | |
|---|---|---|---|---|---|
| <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | |
| 1 | 3113.506 | 4.8994985 | 1634.080 | 3554.460 | 3752.214 |
| 2 | 3260.669 | 2.1678773 | 1595.147 | 3455.778 | 6421.468 |
| 3 | 3416.404 | 2.6359108 | 1559.959 | 3653.636 | 4919.149 |
| 4 | 3407.504 | 2.2857136 | 1570.255 | 3613.218 | 2585.539 |
| 5 | 3489.423 | -1.1614780 | 1582.603 | 3384.890 | 3620.347 |
| 6 | 3303.120 | 0.8239054 | 1648.416 | 3377.272 | 4158.864 |
pred90 %>%
summarize( mean=mean(Ynew), sd=sd(Ynew) )
| mean | sd |
|---|---|
| <dbl> | <dbl> |
| 3515.594 | 1589.728 |
There will be around 3500 bikers!
pred47 <- bike_humidity_model_df %>%
mutate( mu = X.Intercept. + humidity * 47 ) %>%
mutate( Ynew = rnorm( nrow(bike_humidity_model_df), mean=mu, sd=sigma ) )
head( pred47 )
| X.Intercept. | humidity | sigma | mu | Ynew | |
|---|---|---|---|---|---|
| <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | |
| 1 | 3113.506 | 4.8994985 | 1634.080 | 3343.782 | 3033.256 |
| 2 | 3260.669 | 2.1678773 | 1595.147 | 3362.559 | 5396.018 |
| 3 | 3416.404 | 2.6359108 | 1559.959 | 3540.292 | 4599.289 |
| 4 | 3407.504 | 2.2857136 | 1570.255 | 3514.932 | 3470.614 |
| 5 | 3489.423 | -1.1614780 | 1582.603 | 3434.834 | 3320.321 |
| 6 | 3303.120 | 0.8239054 | 1648.416 | 3341.844 | 1736.633 |
pred47 %>%
summarize( mean=mean(Ynew), sd=sd(Ynew) )
| mean | sd |
|---|---|
| <dbl> | <dbl> |
| 3468.608 | 1564.996 |
Also almost 3500 bikers. Our linear model is not very useful here, since the connection between number of riders and humidity is very weak. This makes sense when looking at $\beta_1$ which is centered around zero, resulting in an almost flat line.
pred90$model <- "90%"
pred47$model <- "47%"
pred <- rbind( pred90, pred47 )
ggplot( data=pred, aes(x=Ynew, fill=model) ) + geom_density( alpha=0.4 )
No selectivity whatsoever..
pred47 %>%
summarize( ci_min=quantile(Ynew, 0.1), ci_max=quantile(Ynew, 0.9) )
| ci_min | ci_max |
|---|---|
| <dbl> | <dbl> |
| 1451.187 | 5474.227 |
The number of riders tomorrow will be anywhere between 1451 and 5491 tomorrow with a probability of 80%.
shortcut_prediction90 <- data.frame(
Ynew=posterior_predict(bike_humidity_model, newdata = data.frame(humidity = 90))[,1],
model="90%"
)
shortcut_prediction47 <- data.frame(
Ynew=posterior_predict(bike_humidity_model, newdata = data.frame(humidity = 47))[,1],
model="47%"
)
shortcut_prediction <- rbind( shortcut_prediction90, shortcut_prediction47 )
Distributions:
ggplot( data=shortcut_prediction, aes(x=Ynew, fill=model) ) + geom_density( alpha=0.4 )
A similar result!
80% credible intervals:
shortcut_prediction %>%
group_by( model ) %>%
summarize( ci_min=quantile(Ynew, 0.1), ci_max=quantile(Ynew, 0.9) )
| model | ci_min | ci_max |
|---|---|---|
| <chr> | <dbl> | <dbl> |
| 47% | 1454.840 | 5512.733 |
| 90% | 1459.989 | 5501.419 |
Also here, extremely similar to each other..
The windspeed is given in miles per hour. I am more used to m/s, knots or the Beaufort scale. For simplicity, I use the Beaufort scale to provide rough arguments. Windspeed has probably a correlation with weather conditions and temperature, however this should be modeled separately in a multivariate model.
At average windspeeds, there will probably be an average number of bikers, so I choose the same prior model for the intercept as in the previous exercise ($N(\mu=5000, \sigma=2000)$)). Windspeed probably has only a strong negative impact on the number of bikers when it is strong, however since we deal with a linear model, this effect has to be modeled globally - in general I expect thus a model with a negative slope. Between 0 and 4 Beaufort ($\sim0-15$ mph) it should be pretty ok to bike, however at higher windspeeds $>8$ Beaufort ($\sim 40$ mph) biking is probably no fun anymore, certainly not in storms with $>90$ mph. Assuming the average 5000 bikers at average winds of around 5 mph and no bikers at 90 mph, I get an slope coefficient of
$$\beta_{1,c} \sim \frac{\Delta y}{\Delta x} \sim -\frac{5000 \text{ riders}}{45 \text{ mph}} \sim -111 \text{ riders/mph}$$For extreme values, let us assume that there are 1) still 5000 riders at 90 mph or 2) only 1000 riders at 5 mph, then $\beta_1$ lies somewhere between the extremes $-55$ and $-200$. Tune a normal prior to this end:
options(repr.plot.width=15, repr.plot.height=4)
plot_normal( mean=-111, sd=35 )
I expect the connection between windspeed and number of riders as vague as the connection between humidity and the number of riders and thus assume the same expected value for $\sigma$ as in the previous exercises.
Bayesian Normal regression model structure:
$Y_i = N(\mu_i, \sigma)$, $\mu_i = \beta_0 + \beta_1 X$
$\beta_0 \sim N(\mu=5000, \sigma=2000)$, ; $\beta_{1,c} \sim N(\mu=-111, \sigma=35)$, ; $\sigma \sim \text{Exp}(\lambda=1/2000)$
bike_windspeed_model_prior <- stan_glm(
rides ~ windspeed, data = bikes,
family = gaussian,
prior_intercept = normal(5000, 2000),
prior = normal(-111, 35),
prior_aux = exponential(1/2000),
chains = 5, iter = 4000*2, seed = 84735, prior_PD = TRUE)
SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1). Chain 1: Chain 1: Gradient evaluation took 9e-06 seconds Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds. Chain 1: Adjust your expectations accordingly! Chain 1: Chain 1: Chain 1: Iteration: 1 / 8000 [ 0%] (Warmup) Chain 1: Iteration: 800 / 8000 [ 10%] (Warmup) Chain 1: Iteration: 1600 / 8000 [ 20%] (Warmup) Chain 1: Iteration: 2400 / 8000 [ 30%] (Warmup) Chain 1: Iteration: 3200 / 8000 [ 40%] (Warmup) Chain 1: Iteration: 4000 / 8000 [ 50%] (Warmup) Chain 1: Iteration: 4001 / 8000 [ 50%] (Sampling) Chain 1: Iteration: 4800 / 8000 [ 60%] (Sampling) Chain 1: Iteration: 5600 / 8000 [ 70%] (Sampling) Chain 1: Iteration: 6400 / 8000 [ 80%] (Sampling) Chain 1: Iteration: 7200 / 8000 [ 90%] (Sampling) Chain 1: Iteration: 8000 / 8000 [100%] (Sampling) Chain 1: Chain 1: Elapsed Time: 0.239946 seconds (Warm-up) Chain 1: 0.039221 seconds (Sampling) Chain 1: 0.279167 seconds (Total) Chain 1: SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2). Chain 2: Chain 2: Gradient evaluation took 6e-06 seconds Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds. Chain 2: Adjust your expectations accordingly! Chain 2: Chain 2: Chain 2: Iteration: 1 / 8000 [ 0%] (Warmup) Chain 2: Iteration: 800 / 8000 [ 10%] (Warmup) Chain 2: Iteration: 1600 / 8000 [ 20%] (Warmup) Chain 2: Iteration: 2400 / 8000 [ 30%] (Warmup) Chain 2: Iteration: 3200 / 8000 [ 40%] (Warmup) Chain 2: Iteration: 4000 / 8000 [ 50%] (Warmup) Chain 2: Iteration: 4001 / 8000 [ 50%] (Sampling) Chain 2: Iteration: 4800 / 8000 [ 60%] (Sampling) Chain 2: Iteration: 5600 / 8000 [ 70%] (Sampling) Chain 2: Iteration: 6400 / 8000 [ 80%] (Sampling) Chain 2: Iteration: 7200 / 8000 [ 90%] (Sampling) Chain 2: Iteration: 8000 / 8000 [100%] (Sampling) Chain 2: Chain 2: Elapsed Time: 0.167037 seconds (Warm-up) Chain 2: 0.032271 seconds (Sampling) Chain 2: 0.199308 seconds (Total) Chain 2: SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3). Chain 3: Chain 3: Gradient evaluation took 5e-06 seconds Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds. Chain 3: Adjust your expectations accordingly! Chain 3: Chain 3: Chain 3: Iteration: 1 / 8000 [ 0%] (Warmup) Chain 3: Iteration: 800 / 8000 [ 10%] (Warmup) Chain 3: Iteration: 1600 / 8000 [ 20%] (Warmup) Chain 3: Iteration: 2400 / 8000 [ 30%] (Warmup) Chain 3: Iteration: 3200 / 8000 [ 40%] (Warmup) Chain 3: Iteration: 4000 / 8000 [ 50%] (Warmup) Chain 3: Iteration: 4001 / 8000 [ 50%] (Sampling) Chain 3: Iteration: 4800 / 8000 [ 60%] (Sampling) Chain 3: Iteration: 5600 / 8000 [ 70%] (Sampling) Chain 3: Iteration: 6400 / 8000 [ 80%] (Sampling) Chain 3: Iteration: 7200 / 8000 [ 90%] (Sampling) Chain 3: Iteration: 8000 / 8000 [100%] (Sampling) Chain 3: Chain 3: Elapsed Time: 0.189502 seconds (Warm-up) Chain 3: 0.032291 seconds (Sampling) Chain 3: 0.221793 seconds (Total) Chain 3: SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4). Chain 4: Chain 4: Gradient evaluation took 4e-06 seconds Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds. Chain 4: Adjust your expectations accordingly! Chain 4: Chain 4: Chain 4: Iteration: 1 / 8000 [ 0%] (Warmup) Chain 4: Iteration: 800 / 8000 [ 10%] (Warmup) Chain 4: Iteration: 1600 / 8000 [ 20%] (Warmup) Chain 4: Iteration: 2400 / 8000 [ 30%] (Warmup) Chain 4: Iteration: 3200 / 8000 [ 40%] (Warmup) Chain 4: Iteration: 4000 / 8000 [ 50%] (Warmup) Chain 4: Iteration: 4001 / 8000 [ 50%] (Sampling) Chain 4: Iteration: 4800 / 8000 [ 60%] (Sampling) Chain 4: Iteration: 5600 / 8000 [ 70%] (Sampling) Chain 4: Iteration: 6400 / 8000 [ 80%] (Sampling) Chain 4: Iteration: 7200 / 8000 [ 90%] (Sampling) Chain 4: Iteration: 8000 / 8000 [100%] (Sampling) Chain 4: Chain 4: Elapsed Time: 0.198043 seconds (Warm-up) Chain 4: 0.032879 seconds (Sampling) Chain 4: 0.230922 seconds (Total) Chain 4: SAMPLING FOR MODEL 'continuous' NOW (CHAIN 5). Chain 5: Chain 5: Gradient evaluation took 4e-06 seconds Chain 5: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds. Chain 5: Adjust your expectations accordingly! Chain 5: Chain 5: Chain 5: Iteration: 1 / 8000 [ 0%] (Warmup) Chain 5: Iteration: 800 / 8000 [ 10%] (Warmup) Chain 5: Iteration: 1600 / 8000 [ 20%] (Warmup) Chain 5: Iteration: 2400 / 8000 [ 30%] (Warmup) Chain 5: Iteration: 3200 / 8000 [ 40%] (Warmup) Chain 5: Iteration: 4000 / 8000 [ 50%] (Warmup) Chain 5: Iteration: 4001 / 8000 [ 50%] (Sampling) Chain 5: Iteration: 4800 / 8000 [ 60%] (Sampling) Chain 5: Iteration: 5600 / 8000 [ 70%] (Sampling) Chain 5: Iteration: 6400 / 8000 [ 80%] (Sampling) Chain 5: Iteration: 7200 / 8000 [ 90%] (Sampling) Chain 5: Iteration: 8000 / 8000 [100%] (Sampling) Chain 5: Chain 5: Elapsed Time: 0.319196 seconds (Warm-up) Chain 5: 0.03375 seconds (Sampling) Chain 5: 0.352946 seconds (Total) Chain 5:
Regression lines:
library(tidybayes)
options(repr.plot.width=15, repr.plot.height=10)
bikes %>%
add_fitted_draws(bike_windspeed_model_prior, n = 100) %>%
ggplot(aes(x = windspeed, y = rides)) +
geom_line(aes(y = .value, group = .draw), alpha = 0.1, col = "blue") +
geom_point(data = bikes, size = 0.05)
Warning message: “`fitted_draws` and `add_fitted_draws` are deprecated as their names were confusing. - Use [add_]epred_draws() to get the expectation of the posterior predictive. - Use [add_]linpred_draws() to get the distribution of the linear predictor. - For example, you used [add_]fitted_draws(..., scale = "response"), which means you most likely want [add_]epred_draws(...). NOTE: When updating to the new functions, note that the `model` parameter is now named `object` and the `n` parameter is now named `ndraws`.”
Simulated data:
bikes %>%
add_predicted_draws(bike_windspeed_model_prior, n = 4) %>%
ggplot(aes(x = windspeed, y = rides)) +
geom_point(aes(y = .prediction, group = .draw), size = 0.2) +
facet_wrap(~ .draw)
Warning message:
“
In add_predicted_draws(): The `n` argument is a deprecated alias for `ndraws`.
Use the `ndraws` argument instead.
See help("tidybayes-deprecated").
”
Same as before, intercept varies on a much larger scale than slope. However the intercept is now clearly negative.
ggplot( data=bikes, aes(x=windspeed, y=rides) ) + geom_point() + xlim(0, 35) + ylim(-5000, 17000)
The prior assumptions are a bit too unspecific, the data vary on a much lower scale.
bike_windspeed_model_posterior <- stan_glm(
rides ~ windspeed, data = bikes,
family = gaussian,
prior_intercept = normal(5000, 2000),
prior = normal(-111, 35),
prior_aux = exponential(1/2000),
chains = 5, iter = 4000*2, seed = 84735, prior_PD = FALSE)
SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1). Chain 1: Chain 1: Gradient evaluation took 1.3e-05 seconds Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.13 seconds. Chain 1: Adjust your expectations accordingly! Chain 1: Chain 1: Chain 1: Iteration: 1 / 8000 [ 0%] (Warmup) Chain 1: Iteration: 800 / 8000 [ 10%] (Warmup) Chain 1: Iteration: 1600 / 8000 [ 20%] (Warmup) Chain 1: Iteration: 2400 / 8000 [ 30%] (Warmup) Chain 1: Iteration: 3200 / 8000 [ 40%] (Warmup) Chain 1: Iteration: 4000 / 8000 [ 50%] (Warmup) Chain 1: Iteration: 4001 / 8000 [ 50%] (Sampling) Chain 1: Iteration: 4800 / 8000 [ 60%] (Sampling) Chain 1: Iteration: 5600 / 8000 [ 70%] (Sampling) Chain 1: Iteration: 6400 / 8000 [ 80%] (Sampling) Chain 1: Iteration: 7200 / 8000 [ 90%] (Sampling) Chain 1: Iteration: 8000 / 8000 [100%] (Sampling) Chain 1: Chain 1: Elapsed Time: 0.141302 seconds (Warm-up) Chain 1: 0.143342 seconds (Sampling) Chain 1: 0.284644 seconds (Total) Chain 1: SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2). Chain 2: Chain 2: Gradient evaluation took 9e-06 seconds Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds. Chain 2: Adjust your expectations accordingly! Chain 2: Chain 2: Chain 2: Iteration: 1 / 8000 [ 0%] (Warmup) Chain 2: Iteration: 800 / 8000 [ 10%] (Warmup) Chain 2: Iteration: 1600 / 8000 [ 20%] (Warmup) Chain 2: Iteration: 2400 / 8000 [ 30%] (Warmup) Chain 2: Iteration: 3200 / 8000 [ 40%] (Warmup) Chain 2: Iteration: 4000 / 8000 [ 50%] (Warmup) Chain 2: Iteration: 4001 / 8000 [ 50%] (Sampling) Chain 2: Iteration: 4800 / 8000 [ 60%] (Sampling) Chain 2: Iteration: 5600 / 8000 [ 70%] (Sampling) Chain 2: Iteration: 6400 / 8000 [ 80%] (Sampling) Chain 2: Iteration: 7200 / 8000 [ 90%] (Sampling) Chain 2: Iteration: 8000 / 8000 [100%] (Sampling) Chain 2: Chain 2: Elapsed Time: 0.190084 seconds (Warm-up) Chain 2: 0.143806 seconds (Sampling) Chain 2: 0.33389 seconds (Total) Chain 2: SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3). Chain 3: Chain 3: Gradient evaluation took 8e-06 seconds Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds. Chain 3: Adjust your expectations accordingly! Chain 3: Chain 3: Chain 3: Iteration: 1 / 8000 [ 0%] (Warmup) Chain 3: Iteration: 800 / 8000 [ 10%] (Warmup) Chain 3: Iteration: 1600 / 8000 [ 20%] (Warmup) Chain 3: Iteration: 2400 / 8000 [ 30%] (Warmup) Chain 3: Iteration: 3200 / 8000 [ 40%] (Warmup) Chain 3: Iteration: 4000 / 8000 [ 50%] (Warmup) Chain 3: Iteration: 4001 / 8000 [ 50%] (Sampling) Chain 3: Iteration: 4800 / 8000 [ 60%] (Sampling) Chain 3: Iteration: 5600 / 8000 [ 70%] (Sampling) Chain 3: Iteration: 6400 / 8000 [ 80%] (Sampling) Chain 3: Iteration: 7200 / 8000 [ 90%] (Sampling) Chain 3: Iteration: 8000 / 8000 [100%] (Sampling) Chain 3: Chain 3: Elapsed Time: 0.133841 seconds (Warm-up) Chain 3: 0.142467 seconds (Sampling) Chain 3: 0.276308 seconds (Total) Chain 3: SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4). Chain 4: Chain 4: Gradient evaluation took 7e-06 seconds Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds. Chain 4: Adjust your expectations accordingly! Chain 4: Chain 4: Chain 4: Iteration: 1 / 8000 [ 0%] (Warmup) Chain 4: Iteration: 800 / 8000 [ 10%] (Warmup) Chain 4: Iteration: 1600 / 8000 [ 20%] (Warmup) Chain 4: Iteration: 2400 / 8000 [ 30%] (Warmup) Chain 4: Iteration: 3200 / 8000 [ 40%] (Warmup) Chain 4: Iteration: 4000 / 8000 [ 50%] (Warmup) Chain 4: Iteration: 4001 / 8000 [ 50%] (Sampling) Chain 4: Iteration: 4800 / 8000 [ 60%] (Sampling) Chain 4: Iteration: 5600 / 8000 [ 70%] (Sampling) Chain 4: Iteration: 6400 / 8000 [ 80%] (Sampling) Chain 4: Iteration: 7200 / 8000 [ 90%] (Sampling) Chain 4: Iteration: 8000 / 8000 [100%] (Sampling) Chain 4: Chain 4: Elapsed Time: 0.083435 seconds (Warm-up) Chain 4: 0.145088 seconds (Sampling) Chain 4: 0.228523 seconds (Total) Chain 4: SAMPLING FOR MODEL 'continuous' NOW (CHAIN 5). Chain 5: Chain 5: Gradient evaluation took 7e-06 seconds Chain 5: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds. Chain 5: Adjust your expectations accordingly! Chain 5: Chain 5: Chain 5: Iteration: 1 / 8000 [ 0%] (Warmup) Chain 5: Iteration: 800 / 8000 [ 10%] (Warmup) Chain 5: Iteration: 1600 / 8000 [ 20%] (Warmup) Chain 5: Iteration: 2400 / 8000 [ 30%] (Warmup) Chain 5: Iteration: 3200 / 8000 [ 40%] (Warmup) Chain 5: Iteration: 4000 / 8000 [ 50%] (Warmup) Chain 5: Iteration: 4001 / 8000 [ 50%] (Sampling) Chain 5: Iteration: 4800 / 8000 [ 60%] (Sampling) Chain 5: Iteration: 5600 / 8000 [ 70%] (Sampling) Chain 5: Iteration: 6400 / 8000 [ 80%] (Sampling) Chain 5: Iteration: 7200 / 8000 [ 90%] (Sampling) Chain 5: Iteration: 8000 / 8000 [100%] (Sampling) Chain 5: Chain 5: Elapsed Time: 0.088386 seconds (Warm-up) Chain 5: 0.14425 seconds (Sampling) Chain 5: 0.232636 seconds (Total) Chain 5:
Regression lines (blue) compared with frequentist regression line (red):
library(tidybayes)
options(repr.plot.width=15, repr.plot.height=10)
bikes %>%
add_fitted_draws(bike_windspeed_model_posterior, n = 100) %>%
ggplot(aes(x = windspeed, y = rides)) +
geom_line(aes(y = .value, group = .draw), alpha = 0.1, col = "blue") +
geom_point(data = bikes, size = 0.05) +
ylim(0, 10000) +
geom_smooth( method="lm", col="red" )
Warning message: “`fitted_draws` and `add_fitted_draws` are deprecated as their names were confusing. - Use [add_]epred_draws() to get the expectation of the posterior predictive. - Use [add_]linpred_draws() to get the distribution of the linear predictor. - For example, you used [add_]fitted_draws(..., scale = "response"), which means you most likely want [add_]epred_draws(...). NOTE: When updating to the new functions, note that the `model` parameter is now named `object` and the `n` parameter is now named `ndraws`.” `geom_smooth()` using formula = 'y ~ x'
The variance of the intercept has varied quite a bit!
Evolution of assumed distributions:
bike_windspeed_model_posterior_df = data.frame( bike_windspeed_model_posterior )
bike_windspeed_model_prior_df = data.frame( bike_windspeed_model_prior )
bike_windspeed_model_posterior_df$kind = "posterior"
bike_windspeed_model_prior_df$kind = "prior"
bike_windspeed_model_df = rbind( bike_windspeed_model_prior_df, bike_windspeed_model_posterior_df )
head( bike_windspeed_model_df )
| X.Intercept. | windspeed | sigma | kind | |
|---|---|---|---|---|
| <dbl> | <dbl> | <dbl> | <chr> | |
| 1 | 6842.840 | -103.65511 | 3293.4519 | prior |
| 2 | 5662.131 | -106.60174 | 1785.3671 | prior |
| 3 | 6587.496 | -117.92008 | 1984.4809 | prior |
| 4 | 6862.335 | -128.83726 | 1003.3518 | prior |
| 5 | 6068.603 | -58.25917 | 286.4947 | prior |
| 6 | 5590.360 | -52.30210 | 189.5258 | prior |
Intercept:
options(repr.plot.width=15, repr.plot.height=4)
ggplot( data=bike_windspeed_model_df ) +
geom_density( aes(x=X.Intercept., fill=kind), alpha=0.5 )
Slope:
ggplot( data=bike_windspeed_model_df ) +
geom_density( aes(x=windspeed, fill=kind), alpha=0.5)
Sampling noise:
ggplot( data=bike_windspeed_model_df ) +
geom_density( aes(x=sigma, fill=kind), alpha=0.5)
Our understanding has changed quite a bit. Especially the variance of intercept and sampling noise have reduced drastically. The slope since a bit lower than expected, also with less variance.
head( penguins_bayes )
| species | island | year | bill_length_mm | bill_depth_mm | flipper_length_mm | body_mass_g | above_average_weight | sex |
|---|---|---|---|---|---|---|---|---|
| <fct> | <fct> | <int> | <dbl> | <dbl> | <int> | <int> | <fct> | <fct> |
| Adelie | Torgersen | 2007 | 39.1 | 18.7 | 181 | 3750 | 0 | male |
| Adelie | Torgersen | 2007 | 39.5 | 17.4 | 186 | 3800 | 0 | female |
| Adelie | Torgersen | 2007 | 40.3 | 18.0 | 195 | 3250 | 0 | female |
| Adelie | Torgersen | 2007 | NA | NA | NA | NA | NA | NA |
| Adelie | Torgersen | 2007 | 36.7 | 19.3 | 193 | 3450 | 0 | female |
| Adelie | Torgersen | 2007 | 39.3 | 20.6 | 190 | 3650 | 0 | male |
Prior on bill length:
plot_normal( mean=200, sd=25 )
penguin_model_prior <- stan_glm(
flipper_length_mm ~ bill_length_mm, data = penguins_bayes,
family = gaussian,
prior_intercept = normal(200, 25),
prior = normal(0, 1, autoscale=TRUE),
prior_aux = exponential(1, autoscale=TRUE),
chains = 4, iter = 10000, seed = 84735, prior_PD = TRUE)
SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1). Chain 1: Chain 1: Gradient evaluation took 1e-05 seconds Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds. Chain 1: Adjust your expectations accordingly! Chain 1: Chain 1: Chain 1: Iteration: 1 / 10000 [ 0%] (Warmup) Chain 1: Iteration: 1000 / 10000 [ 10%] (Warmup) Chain 1: Iteration: 2000 / 10000 [ 20%] (Warmup) Chain 1: Iteration: 3000 / 10000 [ 30%] (Warmup) Chain 1: Iteration: 4000 / 10000 [ 40%] (Warmup) Chain 1: Iteration: 5000 / 10000 [ 50%] (Warmup) Chain 1: Iteration: 5001 / 10000 [ 50%] (Sampling) Chain 1: Iteration: 6000 / 10000 [ 60%] (Sampling) Chain 1: Iteration: 7000 / 10000 [ 70%] (Sampling) Chain 1: Iteration: 8000 / 10000 [ 80%] (Sampling) Chain 1: Iteration: 9000 / 10000 [ 90%] (Sampling) Chain 1: Iteration: 10000 / 10000 [100%] (Sampling) Chain 1: Chain 1: Elapsed Time: 0.043339 seconds (Warm-up) Chain 1: 0.03938 seconds (Sampling) Chain 1: 0.082719 seconds (Total) Chain 1: SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2). Chain 2: Chain 2: Gradient evaluation took 7e-06 seconds Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds. Chain 2: Adjust your expectations accordingly! Chain 2: Chain 2: Chain 2: Iteration: 1 / 10000 [ 0%] (Warmup) Chain 2: Iteration: 1000 / 10000 [ 10%] (Warmup) Chain 2: Iteration: 2000 / 10000 [ 20%] (Warmup) Chain 2: Iteration: 3000 / 10000 [ 30%] (Warmup) Chain 2: Iteration: 4000 / 10000 [ 40%] (Warmup) Chain 2: Iteration: 5000 / 10000 [ 50%] (Warmup) Chain 2: Iteration: 5001 / 10000 [ 50%] (Sampling) Chain 2: Iteration: 6000 / 10000 [ 60%] (Sampling) Chain 2: Iteration: 7000 / 10000 [ 70%] (Sampling) Chain 2: Iteration: 8000 / 10000 [ 80%] (Sampling) Chain 2: Iteration: 9000 / 10000 [ 90%] (Sampling) Chain 2: Iteration: 10000 / 10000 [100%] (Sampling) Chain 2: Chain 2: Elapsed Time: 0.040569 seconds (Warm-up) Chain 2: 0.044966 seconds (Sampling) Chain 2: 0.085535 seconds (Total) Chain 2: SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3). Chain 3: Chain 3: Gradient evaluation took 5e-06 seconds Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds. Chain 3: Adjust your expectations accordingly! Chain 3: Chain 3: Chain 3: Iteration: 1 / 10000 [ 0%] (Warmup) Chain 3: Iteration: 1000 / 10000 [ 10%] (Warmup) Chain 3: Iteration: 2000 / 10000 [ 20%] (Warmup) Chain 3: Iteration: 3000 / 10000 [ 30%] (Warmup) Chain 3: Iteration: 4000 / 10000 [ 40%] (Warmup) Chain 3: Iteration: 5000 / 10000 [ 50%] (Warmup) Chain 3: Iteration: 5001 / 10000 [ 50%] (Sampling) Chain 3: Iteration: 6000 / 10000 [ 60%] (Sampling) Chain 3: Iteration: 7000 / 10000 [ 70%] (Sampling) Chain 3: Iteration: 8000 / 10000 [ 80%] (Sampling) Chain 3: Iteration: 9000 / 10000 [ 90%] (Sampling) Chain 3: Iteration: 10000 / 10000 [100%] (Sampling) Chain 3: Chain 3: Elapsed Time: 0.041215 seconds (Warm-up) Chain 3: 0.040103 seconds (Sampling) Chain 3: 0.081318 seconds (Total) Chain 3: SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4). Chain 4: Chain 4: Gradient evaluation took 4e-06 seconds Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds. Chain 4: Adjust your expectations accordingly! Chain 4: Chain 4: Chain 4: Iteration: 1 / 10000 [ 0%] (Warmup) Chain 4: Iteration: 1000 / 10000 [ 10%] (Warmup) Chain 4: Iteration: 2000 / 10000 [ 20%] (Warmup) Chain 4: Iteration: 3000 / 10000 [ 30%] (Warmup) Chain 4: Iteration: 4000 / 10000 [ 40%] (Warmup) Chain 4: Iteration: 5000 / 10000 [ 50%] (Warmup) Chain 4: Iteration: 5001 / 10000 [ 50%] (Sampling) Chain 4: Iteration: 6000 / 10000 [ 60%] (Sampling) Chain 4: Iteration: 7000 / 10000 [ 70%] (Sampling) Chain 4: Iteration: 8000 / 10000 [ 80%] (Sampling) Chain 4: Iteration: 9000 / 10000 [ 90%] (Sampling) Chain 4: Iteration: 10000 / 10000 [100%] (Sampling) Chain 4: Chain 4: Elapsed Time: 0.040952 seconds (Warm-up) Chain 4: 0.038756 seconds (Sampling) Chain 4: 0.079708 seconds (Total) Chain 4:
prior_summary( penguin_model_prior )
Priors for model 'penguin_model_prior'
------
Intercept (after predictors centered)
~ normal(location = 200, scale = 25)
Coefficients
Specified prior:
~ normal(location = 0, scale = 1)
Adjusted prior:
~ normal(location = 0, scale = 2.6)
Auxiliary (sigma)
Specified prior:
~ exponential(rate = 1)
Adjusted prior:
~ exponential(rate = 0.071)
------
See help('prior_summary.stanreg') for more details
Prior model:
$Y_i = N(\mu_i, \sigma)$, $\mu_i = \beta_0 + \beta_1 X$
$\beta_0 \sim N(\mu=200, \sigma=25)$, ; $\beta_{1,c} \sim N(\mu=0, \sigma=2.6)$, ; $\sigma \sim \text{Exp}(\lambda=0.071)$
Regression lines:
penguins_bayes_notna <- penguins_bayes %>%
filter(!is.na(bill_length_mm), !is.na(flipper_length_mm))
library(tidybayes)
options(repr.plot.width=15, repr.plot.height=10)
penguins_bayes_notna %>%
add_fitted_draws(penguin_model_prior, n = 100) %>%
ggplot(aes(x = bill_length_mm, y = flipper_length_mm)) +
geom_line(aes(y = .value, group = .draw), alpha = 0.1, col = "blue") +
geom_point(data = penguins_bayes_notna, size = 0.05)
Warning message: “`fitted_draws` and `add_fitted_draws` are deprecated as their names were confusing. - Use [add_]epred_draws() to get the expectation of the posterior predictive. - Use [add_]linpred_draws() to get the distribution of the linear predictor. - For example, you used [add_]fitted_draws(..., scale = "response"), which means you most likely want [add_]epred_draws(...). NOTE: When updating to the new functions, note that the `model` parameter is now named `object` and the `n` parameter is now named `ndraws`.”
A lot but not anything is possible!
Datasets:
penguins_bayes_notna %>%
add_predicted_draws(penguin_model_prior, n = 4) %>%
ggplot(aes(x = bill_length_mm, y = flipper_length_mm)) +
geom_point(aes(y = .prediction, group = .draw), size = 0.2) +
facet_wrap(~ .draw)
Warning message:
“
In add_predicted_draws(): The `n` argument is a deprecated alias for `ndraws`.
Use the `ndraws` argument instead.
See help("tidybayes-deprecated").
”
options(repr.plot.width=15, repr.plot.height=4)
plot_normal(mean=0, sd=2.6)
x <- seq(0, 60)
plot( x, dexp(x, rate=0.071), type="l" )
Intercepts between 150-250 mm are possible, slopes between -7 and 7 and the noise has a standard deviation of up to 30 mm.
ggplot( data=penguins_bayes, aes(x=bill_length_mm, y=flipper_length_mm) ) + geom_point() + geom_smooth( method="lm")
`geom_smooth()` using formula = 'y ~ x' Warning message: “Removed 2 rows containing non-finite values (`stat_smooth()`).” Warning message: “Removed 2 rows containing missing values (`geom_point()`).”
There is a positive relationship between bill and flipper length with a slope around 1.5, the standard deviation of noise is $<10$ mm.
Yes, the residuals look more or less normally distributed and the variance more or less constant.
options(repr.plot.width=15, repr.plot.height=4)
ggplot( penguins_bayes %>% filter( bill_length_mm>35, bill_length_mm<45 ), aes(x=flipper_length_mm) ) +
geom_histogram()
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
penguin_model_posterior <- stan_glm(
flipper_length_mm ~ bill_length_mm, data = penguins_bayes,
family = gaussian,
prior_intercept = normal(200, 25),
prior = normal(0, 1, autoscale=TRUE),
prior_aux = exponential(1, autoscale=TRUE),
chains = 4, iter = 10000, seed = 84735, prior_PD = FALSE)
SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1). Chain 1: Chain 1: Gradient evaluation took 1.4e-05 seconds Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds. Chain 1: Adjust your expectations accordingly! Chain 1: Chain 1: Chain 1: Iteration: 1 / 10000 [ 0%] (Warmup) Chain 1: Iteration: 1000 / 10000 [ 10%] (Warmup) Chain 1: Iteration: 2000 / 10000 [ 20%] (Warmup) Chain 1: Iteration: 3000 / 10000 [ 30%] (Warmup) Chain 1: Iteration: 4000 / 10000 [ 40%] (Warmup) Chain 1: Iteration: 5000 / 10000 [ 50%] (Warmup) Chain 1: Iteration: 5001 / 10000 [ 50%] (Sampling) Chain 1: Iteration: 6000 / 10000 [ 60%] (Sampling) Chain 1: Iteration: 7000 / 10000 [ 70%] (Sampling) Chain 1: Iteration: 8000 / 10000 [ 80%] (Sampling) Chain 1: Iteration: 9000 / 10000 [ 90%] (Sampling) Chain 1: Iteration: 10000 / 10000 [100%] (Sampling) Chain 1: Chain 1: Elapsed Time: 0.059852 seconds (Warm-up) Chain 1: 0.139061 seconds (Sampling) Chain 1: 0.198913 seconds (Total) Chain 1: SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2). Chain 2: Chain 2: Gradient evaluation took 9e-06 seconds Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds. Chain 2: Adjust your expectations accordingly! Chain 2: Chain 2: Chain 2: Iteration: 1 / 10000 [ 0%] (Warmup) Chain 2: Iteration: 1000 / 10000 [ 10%] (Warmup) Chain 2: Iteration: 2000 / 10000 [ 20%] (Warmup) Chain 2: Iteration: 3000 / 10000 [ 30%] (Warmup) Chain 2: Iteration: 4000 / 10000 [ 40%] (Warmup) Chain 2: Iteration: 5000 / 10000 [ 50%] (Warmup) Chain 2: Iteration: 5001 / 10000 [ 50%] (Sampling) Chain 2: Iteration: 6000 / 10000 [ 60%] (Sampling) Chain 2: Iteration: 7000 / 10000 [ 70%] (Sampling) Chain 2: Iteration: 8000 / 10000 [ 80%] (Sampling) Chain 2: Iteration: 9000 / 10000 [ 90%] (Sampling) Chain 2: Iteration: 10000 / 10000 [100%] (Sampling) Chain 2: Chain 2: Elapsed Time: 0.055711 seconds (Warm-up) Chain 2: 0.14411 seconds (Sampling) Chain 2: 0.199821 seconds (Total) Chain 2: SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3). Chain 3: Chain 3: Gradient evaluation took 6e-06 seconds Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds. Chain 3: Adjust your expectations accordingly! Chain 3: Chain 3: Chain 3: Iteration: 1 / 10000 [ 0%] (Warmup) Chain 3: Iteration: 1000 / 10000 [ 10%] (Warmup) Chain 3: Iteration: 2000 / 10000 [ 20%] (Warmup) Chain 3: Iteration: 3000 / 10000 [ 30%] (Warmup) Chain 3: Iteration: 4000 / 10000 [ 40%] (Warmup) Chain 3: Iteration: 5000 / 10000 [ 50%] (Warmup) Chain 3: Iteration: 5001 / 10000 [ 50%] (Sampling) Chain 3: Iteration: 6000 / 10000 [ 60%] (Sampling) Chain 3: Iteration: 7000 / 10000 [ 70%] (Sampling) Chain 3: Iteration: 8000 / 10000 [ 80%] (Sampling) Chain 3: Iteration: 9000 / 10000 [ 90%] (Sampling) Chain 3: Iteration: 10000 / 10000 [100%] (Sampling) Chain 3: Chain 3: Elapsed Time: 0.053597 seconds (Warm-up) Chain 3: 0.142492 seconds (Sampling) Chain 3: 0.196089 seconds (Total) Chain 3: SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4). Chain 4: Chain 4: Gradient evaluation took 7e-06 seconds Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds. Chain 4: Adjust your expectations accordingly! Chain 4: Chain 4: Chain 4: Iteration: 1 / 10000 [ 0%] (Warmup) Chain 4: Iteration: 1000 / 10000 [ 10%] (Warmup) Chain 4: Iteration: 2000 / 10000 [ 20%] (Warmup) Chain 4: Iteration: 3000 / 10000 [ 30%] (Warmup) Chain 4: Iteration: 4000 / 10000 [ 40%] (Warmup) Chain 4: Iteration: 5000 / 10000 [ 50%] (Warmup) Chain 4: Iteration: 5001 / 10000 [ 50%] (Sampling) Chain 4: Iteration: 6000 / 10000 [ 60%] (Sampling) Chain 4: Iteration: 7000 / 10000 [ 70%] (Sampling) Chain 4: Iteration: 8000 / 10000 [ 80%] (Sampling) Chain 4: Iteration: 9000 / 10000 [ 90%] (Sampling) Chain 4: Iteration: 10000 / 10000 [100%] (Sampling) Chain 4: Chain 4: Elapsed Time: 0.057119 seconds (Warm-up) Chain 4: 0.143422 seconds (Sampling) Chain 4: 0.200541 seconds (Total) Chain 4:
library(tidybayes)
options(repr.plot.width=15, repr.plot.height=10)
penguins_bayes_notna %>%
add_fitted_draws(penguin_model_posterior, n = 100) %>%
ggplot(aes(x = bill_length_mm, y = flipper_length_mm)) +
geom_line(aes(y = .value, group = .draw), alpha = 0.1, col = "blue") +
geom_point(data = penguins_bayes_notna, size = 0.05)
Warning message: “`fitted_draws` and `add_fitted_draws` are deprecated as their names were confusing. - Use [add_]epred_draws() to get the expectation of the posterior predictive. - Use [add_]linpred_draws() to get the distribution of the linear predictor. - For example, you used [add_]fitted_draws(..., scale = "response"), which means you most likely want [add_]epred_draws(...). NOTE: When updating to the new functions, note that the `model` parameter is now named `object` and the `n` parameter is now named `ndraws`.”
The data have narrowed down the distribution for the slope quite a bit.
tidy(penguin_model_posterior, effects = c("fixed", "aux"),
conf.int = TRUE, conf.level = 0.9)
| term | estimate | std.error | conf.low | conf.high |
|---|---|---|---|---|
| <chr> | <dbl> | <dbl> | <dbl> | <dbl> |
| (Intercept) | 126.782719 | 4.7034754 | 119.137746 | 134.400136 |
| bill_length_mm | 1.687882 | 0.1066332 | 1.516402 | 1.860831 |
| sigma | 10.637865 | 0.4102912 | 10.005861 | 11.342029 |
| mean_PPD | 200.923677 | 0.8190803 | 199.570673 | 202.259487 |
With 90% plausibility, the slope $\beta_1$ lies in the interval [1.52,1.86].
Yes, we have ample evidence. Penguins with longer bills tend to have longer flippers as long as $\beta_1 > 0$. In the posterior simulation, all coefficients are larger than zero, so we currently believe that there is a positive relationship with a 100% plausibility.
penguin_model_posterior_df = data.frame(penguin_model_posterior)
penguin_model_posterior_df %>%
summarize( fraction=mean(bill_length_mm>0) )
| fraction |
|---|
| <dbl> |
| 1 |
Simulation:
sim <- penguin_model_posterior_df %>%
mutate( mu=X.Intercept. + 51*bill_length_mm ) %>%
mutate( Ynew=rnorm(nrow(penguin_model_posterior_df), mean=mu, sd=sigma) )
head( sim )
| X.Intercept. | bill_length_mm | sigma | mu | Ynew | |
|---|---|---|---|---|---|
| <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | |
| 1 | 123.8772 | 1.759952 | 11.17994 | 213.6348 | 219.7600 |
| 2 | 124.4226 | 1.729037 | 10.06938 | 212.6035 | 229.1960 |
| 3 | 129.5158 | 1.639111 | 11.28116 | 213.1105 | 224.1757 |
| 4 | 126.0923 | 1.717777 | 10.46458 | 213.6989 | 216.4598 |
| 5 | 128.1085 | 1.662646 | 10.88810 | 212.9035 | 207.0505 |
| 6 | 130.6763 | 1.621520 | 10.30822 | 213.3738 | 215.3477 |
Posterior model for typical flipper length:
sim %>%
summarize( lower=quantile(mu, 0.05), upper=quantile(mu,0.95) )
| lower | upper |
|---|---|
| <dbl> | <dbl> |
| 211.3178 | 214.4123 |
The average flipper length for all penguins with bill length of 51 mm lies with a probability of 90% between 211.3 and 214.4.
Posterior predictive model:
sim %>%
summarize( lower=quantile(Ynew, 0.05), upper=quantile(Ynew,0.95) )
| lower | upper |
|---|---|
| <dbl> | <dbl> |
| 195.5244 | 230.2396 |
Pablo's flipper length lies between 195.3 and 230.4 mm with a probability of 90%. The contribution of sampling variability is quite large here.
Average flipper posterior length in red, posterior predictive model in blue.
options(repr.plot.width=15, repr.plot.height=4)
ggplot( sim ) +
geom_density( aes(x=mu), fill="red", alpha=0.5 ) +
geom_density( aes(x=Ynew), fill="blue", alpha=0.5 )
sim %>%
summarize( lower=quantile(Ynew, 0.1), upper=quantile(Ynew,0.9) )
| lower | upper |
|---|---|
| <dbl> | <dbl> |
| 199.2137 | 226.537 |
More narrow, since sampling variability has not to be accounted for. Here's the proof:
sim %>%
summarize( lower=quantile(mu, 0.1), upper=quantile(mu,0.9) )
| lower | upper |
|---|---|
| <dbl> | <dbl> |
| 211.6633 | 214.0649 |
shortcut_prediction <- data.frame(
Ynew=posterior_predict(penguin_model_posterior, newdata = data.frame(bill_length_mm = 51))[,1]
)
shortcut_prediction %>%
summarize( lower=quantile(Ynew, 0.1), upper=quantile(Ynew,0.9) )
| lower | upper |
|---|---|
| <dbl> | <dbl> |
| 199.1529 | 226.6365 |
Practically similar!
options(repr.plot.width=15, repr.plot.height=4)
ggplot( shortcut_prediction ) +
geom_density( aes(x=Ynew), fill="blue", alpha=0.5 )
plot_normal( mean=0.01, sd=0.002 )
qnorm( p=c(0.05, 0.95) , mean=0.01, sd=0.002 )
Researchers are quite certain, that the slope of the relationship is somewhere between 0.005 and 0.015. In other words, they believe to 90% that the slope is somewhere between 0.0067 and 0.013.
ggplot( data=penguins_bayes, aes(x=body_mass_g, y=flipper_length_mm) ) + geom_point() + geom_smooth( method="lm")
`geom_smooth()` using formula = 'y ~ x' Warning message: “Removed 2 rows containing non-finite values (`stat_smooth()`).” Warning message: “Removed 2 rows containing missing values (`geom_point()`).”
Observed relationship: Slope of around 0.016, larger than expected.
In the posterior model using bill_length, we found $\sigma \sim 10$. It appears that $\sigma$ here is slightly lower, but not much more than a factor of two.
Estimate of around 7:
penguins_bayes %>%
filter( body_mass_g>3700, body_mass_g<4000 ) %>%
summarise( mean=mean(flipper_length_mm), sd=sd(flipper_length_mm) )
| mean | sd |
|---|---|
| <dbl> | <dbl> |
| 193.1042 | 7.131708 |
penguin_model2_posterior <- stan_glm(
flipper_length_mm ~ body_mass_g, data = penguins_bayes,
family = gaussian,
prior_intercept = normal(200, 25, autoscale=TRUE),
prior = normal(0.01, 0.002),
prior_aux = exponential(1, autoscale=TRUE),
chains = 4, iter = 10000, seed = 84735, prior_PD = FALSE)
SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1). Chain 1: Chain 1: Gradient evaluation took 1.3e-05 seconds Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.13 seconds. Chain 1: Adjust your expectations accordingly! Chain 1: Chain 1: Chain 1: Iteration: 1 / 10000 [ 0%] (Warmup) Chain 1: Iteration: 1000 / 10000 [ 10%] (Warmup) Chain 1: Iteration: 2000 / 10000 [ 20%] (Warmup) Chain 1: Iteration: 3000 / 10000 [ 30%] (Warmup) Chain 1: Iteration: 4000 / 10000 [ 40%] (Warmup) Chain 1: Iteration: 5000 / 10000 [ 50%] (Warmup) Chain 1: Iteration: 5001 / 10000 [ 50%] (Sampling) Chain 1: Iteration: 6000 / 10000 [ 60%] (Sampling) Chain 1: Iteration: 7000 / 10000 [ 70%] (Sampling) Chain 1: Iteration: 8000 / 10000 [ 80%] (Sampling) Chain 1: Iteration: 9000 / 10000 [ 90%] (Sampling) Chain 1: Iteration: 10000 / 10000 [100%] (Sampling) Chain 1: Chain 1: Elapsed Time: 0.058577 seconds (Warm-up) Chain 1: 0.143265 seconds (Sampling) Chain 1: 0.201842 seconds (Total) Chain 1: SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2). Chain 2: Chain 2: Gradient evaluation took 8e-06 seconds Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds. Chain 2: Adjust your expectations accordingly! Chain 2: Chain 2: Chain 2: Iteration: 1 / 10000 [ 0%] (Warmup) Chain 2: Iteration: 1000 / 10000 [ 10%] (Warmup) Chain 2: Iteration: 2000 / 10000 [ 20%] (Warmup) Chain 2: Iteration: 3000 / 10000 [ 30%] (Warmup) Chain 2: Iteration: 4000 / 10000 [ 40%] (Warmup) Chain 2: Iteration: 5000 / 10000 [ 50%] (Warmup) Chain 2: Iteration: 5001 / 10000 [ 50%] (Sampling) Chain 2: Iteration: 6000 / 10000 [ 60%] (Sampling) Chain 2: Iteration: 7000 / 10000 [ 70%] (Sampling) Chain 2: Iteration: 8000 / 10000 [ 80%] (Sampling) Chain 2: Iteration: 9000 / 10000 [ 90%] (Sampling) Chain 2: Iteration: 10000 / 10000 [100%] (Sampling) Chain 2: Chain 2: Elapsed Time: 0.068671 seconds (Warm-up) Chain 2: 0.142196 seconds (Sampling) Chain 2: 0.210867 seconds (Total) Chain 2: SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3). Chain 3: Chain 3: Gradient evaluation took 8e-06 seconds Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds. Chain 3: Adjust your expectations accordingly! Chain 3: Chain 3: Chain 3: Iteration: 1 / 10000 [ 0%] (Warmup) Chain 3: Iteration: 1000 / 10000 [ 10%] (Warmup) Chain 3: Iteration: 2000 / 10000 [ 20%] (Warmup) Chain 3: Iteration: 3000 / 10000 [ 30%] (Warmup) Chain 3: Iteration: 4000 / 10000 [ 40%] (Warmup) Chain 3: Iteration: 5000 / 10000 [ 50%] (Warmup) Chain 3: Iteration: 5001 / 10000 [ 50%] (Sampling) Chain 3: Iteration: 6000 / 10000 [ 60%] (Sampling) Chain 3: Iteration: 7000 / 10000 [ 70%] (Sampling) Chain 3: Iteration: 8000 / 10000 [ 80%] (Sampling) Chain 3: Iteration: 9000 / 10000 [ 90%] (Sampling) Chain 3: Iteration: 10000 / 10000 [100%] (Sampling) Chain 3: Chain 3: Elapsed Time: 0.05902 seconds (Warm-up) Chain 3: 0.141205 seconds (Sampling) Chain 3: 0.200225 seconds (Total) Chain 3: SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4). Chain 4: Chain 4: Gradient evaluation took 5e-06 seconds Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds. Chain 4: Adjust your expectations accordingly! Chain 4: Chain 4: Chain 4: Iteration: 1 / 10000 [ 0%] (Warmup) Chain 4: Iteration: 1000 / 10000 [ 10%] (Warmup) Chain 4: Iteration: 2000 / 10000 [ 20%] (Warmup) Chain 4: Iteration: 3000 / 10000 [ 30%] (Warmup) Chain 4: Iteration: 4000 / 10000 [ 40%] (Warmup) Chain 4: Iteration: 5000 / 10000 [ 50%] (Warmup) Chain 4: Iteration: 5001 / 10000 [ 50%] (Sampling) Chain 4: Iteration: 6000 / 10000 [ 60%] (Sampling) Chain 4: Iteration: 7000 / 10000 [ 70%] (Sampling) Chain 4: Iteration: 8000 / 10000 [ 80%] (Sampling) Chain 4: Iteration: 9000 / 10000 [ 90%] (Sampling) Chain 4: Iteration: 10000 / 10000 [100%] (Sampling) Chain 4: Chain 4: Elapsed Time: 0.055587 seconds (Warm-up) Chain 4: 0.142819 seconds (Sampling) Chain 4: 0.198406 seconds (Total) Chain 4:
penguin_model2_posterior_df = data.frame( penguin_model2_posterior )
ggplot( data=penguin_model2_posterior_df, aes( x=body_mass_g ) ) +
geom_density( fill="blue", alpha=0.5 ) +
stat_function(fun = dnorm, args = list(mean = 0.01, sd = 0.002), color="red") + ylab("") +
scale_y_continuous(breaks = NULL) +
xlim(0.005, 0.02)
quantile( penguin_model2_posterior_df$body_mass_g, c(0.05,0.95) )
Variance has decreased and the mean increased significantly. The researchers can successfully update their prior. Before they assumed with a probability of 90% that $\beta_1$ is somewhere between 0.0067 and 0.0134, now this interval has narrowed to somewhere between 0.0142 and 0.0157.